function [mergers,subs,maxs,nA,nF,Omega]=equilibrium(sigmaS,sigmaT,N,T,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nSim)

mu = 1/T;

%Values
L_sp_end=1;
L_team_end=1/2;
F_end = 0;
L_sp=zeros(N/2+1,nScore);
F_sp=zeros(N/2+1,nScore,2);
L_team=zeros(N/2+1,nScore,1);
F_team=F_sp;
psi=1-lambda1-lambda2+2*lambda2*[0:1:N/2]'/N;
Psub_f_team=F_sp;
Psub_f_sp=F_sp;
Psub_l_team=L_sp;
Psub_l_sp=L_sp;
Pteam=F_sp;
Exit_Flag=L_sp;

L_team(:,end)=L_team_end;
L_sp(:,end)=L_sp_end;


%%

K = @(x,sigma) x.^sigma;
EC_Z = @(z,sigma) z*sigma/(1+sigma); %E[c|c<z] (truncated expectation)
fun_p = @(q,V1,V2,sigma) K(max(q*(V1-V2),0),sigma);
fun_exp_max = @(p,q,V1,V2,sigma) p*(q*V1+(1-q)*V2+EC_Z(q*(V1-V2),sigma))+(1-p)*V2; 
options = optimset('Display','off','MaxIter',50000,'MaxFunEvals',50000,...
            'TolX', 1e-9, 'TolFun', 1e-7);
for i=0:(N/2)
    aux_Teams = N/2 - i;
    iaux_Teams = aux_Teams+1;
    aux_sp=N - 2*aux_Teams;
    for j=2:nScore
        aux_Score = nScore - j + 1;
        
        %x0 = [F_sp(aux_Teams,aux_Score+1,1+0);F_sp(aux_Teams,aux_Score+1,1+1);F_team(aux_Teams,aux_Score+1,1+0);F_team(aux_Teams,aux_Score+1,1+1);L_sp(aux_Teams,aux_Score+1); L_team(aux_Teams,aux_Score+1)];
        x0=zeros(6,1);
        [xstar,error_star,exit_flag] = fsolve(@(x) equilibrium_mapping(x,mu,sigmaS,sigmaT,N,lambda1,lambda2,q_team,q_sp,L_sp_end,L_team_end,F_end,L_sp,F_sp,L_team,F_team,aux_Teams,aux_Score,aux_sp),x0,options);
       
        Exit_Flag(iaux_Teams,aux_Score)=exit_flag;
        xstar=max(xstar,0);

%    x0=zeros(6,1);
%    error=1;
%     while error>1e-6
%         x1=equilibrium_mapping(x0,mu,sigmaS,sigmaT,N,lambda1,lambda2,q_team,q_sp,L_sp_end,L_team_end,F_end,L_sp,F_sp,L_team,F_team,aux_Teams,aux_Score,aux_sp);
%         error=(x1-x0)'*(x1-x0);
%         x0=x1;
%     end

        F_sp_0=xstar(1);
        F_sp_1=xstar(2);
        F_team_0=xstar(3);
        F_team_1=xstar(4);
        L_sp_0=xstar(5);
        L_team_1=xstar(6);
        
        p_sub_f_sp_0 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1),F_sp(iaux_Teams,aux_Score,1+0),sigmaS);
        p_sub_f_sp_1 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1),F_sp(iaux_Teams,aux_Score,1+1),sigmaS);
        p_sub_f_team_0 = fun_p(q_team(aux_Score),L_team(iaux_Teams,aux_Score+1),F_team(iaux_Teams,aux_Score,1+0),sigmaS);
        p_sub_f_team_1 = fun_p(q_team(aux_Score),L_team(iaux_Teams,aux_Score+1),F_team(iaux_Teams,aux_Score,1+1),sigmaS);        
        p_sub_l_sp_0 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1),L_sp(iaux_Teams,aux_Score),sigmaS);
        p_sub_l_team_1 = fun_p(q_team(aux_Score),L_team(iaux_Teams,aux_Score+1),L_team(iaux_Teams,aux_Score),sigmaS);
        if aux_Teams<N/2
            p_team_f_sp_0 = fun_p(1,F_sp(iaux_Teams+1,aux_Score),F_sp(iaux_Teams,aux_Score,1+0),sigmaT);
            p_team_f_sp_1 = fun_p(1,F_sp(iaux_Teams+1,aux_Score),F_sp(iaux_Teams,aux_Score,1+1),sigmaT);
        else
            p_team_f_sp_0=0;
            p_team_f_sp_1=0;
        end

        Psub_f_team(iaux_Teams,aux_Score,1)=p_sub_f_team_0;
        Psub_f_team(iaux_Teams,aux_Score,2)=p_sub_f_team_1;
        Psub_f_sp(iaux_Teams,aux_Score,1)=p_sub_f_sp_0;
        Psub_f_sp(iaux_Teams,aux_Score,2)=p_sub_f_sp_1;
        Psub_l_team(iaux_Teams,aux_Score)=p_sub_l_team_1;
        Psub_l_sp(iaux_Teams,aux_Score)=p_sub_l_sp_0;
        Pteam(iaux_Teams,aux_Score,1)=p_team_f_sp_0;
        Pteam(iaux_Teams,aux_Score,2)=p_team_f_sp_1;
        
        F_sp(iaux_Teams,aux_Score,1)=F_sp_0;
        F_sp(iaux_Teams,aux_Score,2)=F_sp_1;
        F_team(iaux_Teams,aux_Score,1)=F_team_0;
        F_team(iaux_Teams,aux_Score,2)=F_team_1;  
        L_sp(iaux_Teams,aux_Score)=L_sp_0;
        L_team(iaux_Teams,aux_Score)=L_team_1;
        
          end
end


%%
Out=zeros(nSim,6);
parfor k=1:nSim
    state=ones(T,2);
    subs=zeros(T,1);
    mergers=zeros(T,1);
    for t=1:T
        rng(t+T*k)
        subs(t,1)=unifrnd(0,1)<=lambda1;
        nA = 2*matN(state(t,1),1);
        nSP = N - 2*sum(matN(state(t,1),:));
        if subs(t,1)==1
            aux1 = unifrnd(0,1);
            aux2 = unifrnd(0,1);
            state(t+1,1) = state(t,1);
            if aux1<nA/N
                if aux2<q_team(state(t,2),1)
                    state(t+1,2) = state(t,2) + 1;
                else
                    state(t+1,2) = state(t,2);
                end
            end
            if aux1>nA/N && aux1<(nSP+nA)/N
                if aux2<q_sp(state(t,2),1)
                    state(t+1,2) = state(t,2) + 1;
                else
                    state(t+1,2) = state(t,2);
                end
            end
            if aux1>(nSP+nA)/N
                state(t+1,2) = state(t,2);
                subs(t,1)=0;
            end
        else
            aux1 = unifrnd(0,1);
            aux2 = unifrnd(0,1);
            aux3 = unifrnd(0,1);
            state(t+1,2) = state(t,2);
            if aux1<=nSP/N
                if aux2<=P(state(t,1),state(t,2))
                    mergers(t,1)=1;
                    if aux3<gamma
                        state(t+1,1) = tran_n(state(t,1),1);
                    else
                        state(t+1,1) = tran_n(state(t,1),2);
                    end
                else
                    state(t+1,1) = state(t,1);
                end
            else
                state(t+1,1) = state(t,1);
            end
        end
    end
    state=state(1:T,:);
    Out(k,:)=[sum(subs),Svalues(state(T,2)),sum(mergers),matN(state(T,1),:),(sum(mergers)-mergers_data).^2];
end
Out=mean(Out,1);
subs=Out(1);
maxs=Out(2);
mergers=Out(3);
nA=Out(4);
nF=Out(5);
Omega=Out(6);
end
